library(raster)

library(plyr)
library(dplyr)
library(reshape)

library(ggrepel)
library(ggplot2)
library(ggpmisc)
library(ggpubr)

library(viridisLite)
library(colorspace)
library(RColorBrewer)

library(sf)
library(stringr)

library(kableExtra)

memory.limit(size = 160000)
## [1] Inf

#display.brewer.pal(n = 11, name ="RdYlGn")
colorvec <- brewer.pal(n = 11, name ="RdYlGn")

col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "NA" = "#D5D8DC")

col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

1 Data

Liens utiles:

# era5vA <- read.table("era5_AI.txt") %>%
#   mutate(tas = tas-273.15, pr = mpr * 60 * 60 * 24 * 365, rsds = ssrd / (60*60*24) * 1e-6 * 60*60*24) %>% # conversion from J/m2 to J/sm2 to MJ/m2/day  
#   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","rsds","ea","ET0","AI","cat.AI"))
# 
# 
# write.table(era5vA, "era5vA.txt")

era5vA <- read.table("ERA5/era5_AI.txt") %>% 
  mutate(z = z, tas = t2m - 273.15, sfcWind = wind, model = "historical", source = "ERA5") %>%  #pr, ET0 already in mm/year
  select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","ea","z","ET0","AI","cat.AI")) %>%
  filter(lm == 1)

write.table(era5vA, "era5vA.txt")

wdcvA <- read.table("wdc_AI.txt") %>% mutate(z = z, ea = vapr, ET0 = ET0) %>% # srad in MJ/m2/day, pr and ETO already in mm/year
 select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI","cat.AI")) %>%
  filter(lm == 1)


write.table(wdcvA, "wdcvA.txt")


cmip6 <- read.table("cmip6_AI.txt")

cmip6vA <- cmip6 %>% subset(period == "1970_2000") %>%
   group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, z) %>% 
   dplyr::summarise(tas = mean(tas, na.rm = T) - 273.15, pr = mean(spr, na.rm = T), # pr in mm/year
                   sfcWind = mean(sfcWind, na.rm = T), Rn = mean(Rn, na.rm = T),
                   ea = mean(ea, na.rm = T), ET0 = mean(sET0, na.rm = T), AI = mean(AI, na.rm = T)) %>% # ET0 in mm/year
  ungroup() %>%  
  mutate(source = "CMIP6")%>%
   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI"))


breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

cmip6vA$cat.AI <- cmip6vA$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6vA$cat.AI[which(cmip6vA$ET0 < 400)] <- "Cold"

write.table(cmip6vA, "cmip6vA.txt")
era5vA <- read.table("figures_article/era5vA.txt")
wdcvA <- read.table("figures_article/wdcvA.txt")
cmip6vA <- read.table("figures_article/cmip6vA.txt")

2 Regions and Continents


colors_continents <- c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
                  "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
                  "#aec7e8", "#ffbb78", "#98df8a", "#ff9896")

ggplot(era5vA)+geom_raster(aes(x= lon, y = lat, fill = Continent))+
  scale_fill_discrete(type = colors_continents)+
  theme_void()

“Continents” to remove are: * Arctic (too little gridcells) * Atlantic (ocean) * Indian (ocean) * Pacific (ocean) * Southern (ocean)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
plot(land_mask)

land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land


elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif") 
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))


ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326")) 

ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))

ipcc_sf <- st_as_sf(ipcc_regions)

ipcc_sf$Name2 <- ipcc_sf$Name %>% str_wrap(width = 10)

ggplot(data = filter(ipcc_sf, Type %in% c("Land","Land-Ocean")))+
  geom_sf(aes(fill=Continent))+
  scale_fill_manual(values = colorvec)+
 geom_sf_text(aes(label = Name), size = 3, position = position_jitter(height = 3, width = 3, seed = 1))+
  #geom_sf_text(aes(label = Acronym), size = 2)+
  theme_void()+
  theme(legend.position = "bottom")

#ggsave("ipcc_regions.png", width = 12, height = 8, units = "cm", dpi = "retina")

3 Aridity Index FAO

3.1 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000 with country borders


calib <- rbind(cmip6vA, era5vA, wdcvA)

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"

calib$cat.AI <- factor(calib$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top"),

nrows = 4, ncol = 1

)

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggsave("worldclim_AIref.png", width = 12, height = 8, dpi = "retina")

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggsave("era5_AIref.png", width = 12, height = 8, dpi = "retina")

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggsave("cmip6_AIref.png", width = 12, height = 8, dpi = "retina")
wdcfordiff <- calib %>% subset(source == "Worldclim") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.wdc")) 
era5fordiff <- calib %>% subset(source == "ERA5") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.era5")) 
cmip6fordiff <- calib %>% subset(source == "CMIP6") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmip6")) 

diffcatAI <- wdcfordiff %>% merge(era5fordiff, by = c("lon","lat")) %>% merge(cmip6fordiff, by = c("lon","lat"))

diffcatAI$era5_wdc <- with(diffcatAI, 
          ifelse(cat.AI.era5 == "Arid" & cat.AI.wdc %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
          ifelse(cat.AI.era5 == "Arid" & cat.AI.wdc == "Hyperarid", "dryer",
          ifelse(cat.AI.era5 == "Semi-arid" & cat.AI.wdc %in% c("Cold","Humid","Dry subhumid"), "wetter",
          ifelse(cat.AI.era5 == "Semi-arid" & cat.AI.wdc %in% c("Hyperarid","Arid"), "dryer",
          ifelse(cat.AI.era5 == "Dry subhumid" & cat.AI.wdc %in% c("Cold","Humid"), "wetter",
          ifelse(cat.AI.era5 == "Dry subhumid" & cat.AI.wdc %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Humid" & cat.AI.wdc == "Cold","wetter",
          ifelse(cat.AI.era5 == "Humid" & cat.AI.wdc %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Cold" & cat.AI.wdc %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Hyperarid" & cat.AI.wdc %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
          ifelse(cat.AI.era5 == cat.AI.wdc, "same",
          NA))))))))))))

diffcatAI$era5_cmip6 <- with(diffcatAI, 
          ifelse(cat.AI.era5 == "Arid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
          ifelse(cat.AI.era5 == "Arid" & cat.AI.cmip6 == "Hyperarid", "dryer",
          ifelse(cat.AI.era5 == "Semi-arid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid"), "wetter",
          ifelse(cat.AI.era5 == "Semi-arid" & cat.AI.cmip6 %in% c("Hyperarid","Arid"), "dryer",
          ifelse(cat.AI.era5 == "Dry subhumid" & cat.AI.cmip6 %in% c("Cold","Humid"), "wetter",
          ifelse(cat.AI.era5 == "Dry subhumid" & cat.AI.cmip6 %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Humid" & cat.AI.cmip6 == "Cold","wetter",
          ifelse(cat.AI.era5 == "Humid" & cat.AI.cmip6 %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Cold" & cat.AI.cmip6 %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Hyperarid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
          ifelse(cat.AI.era5 == cat.AI.cmip6, "same",
          NA))))))))))))

diffcatAI$wdc_cmip6 <- with(diffcatAI, 
          ifelse(cat.AI.wdc == "Arid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
          ifelse(cat.AI.wdc == "Arid" & cat.AI.cmip6 == "Hyperarid", "dryer",
          ifelse(cat.AI.wdc == "Semi-arid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid"), "wetter",
          ifelse(cat.AI.wdc == "Semi-arid" & cat.AI.cmip6 %in% c("Hyperarid","Arid"), "dryer",
          ifelse(cat.AI.wdc == "Dry subhumid" & cat.AI.cmip6 %in% c("Cold","Humid"), "wetter",
          ifelse(cat.AI.wdc == "Dry subhumid" & cat.AI.cmip6 %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.wdc == "Humid" & cat.AI.cmip6 == "Cold","wetter",
          ifelse(cat.AI.wdc == "Humid" & cat.AI.cmip6 %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.wdc == "Cold" & cat.AI.cmip6 %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.wdc == "Hyperarid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
          ifelse(cat.AI.wdc == cat.AI.cmip6, "same",
          NA))))))))))))

ggplot(filter(diffcatAI, cat.AI.era5 != cat.AI.wdc)) + geom_raster(aes(x=lon, y = lat,  fill = era5_wdc))+
    scale_fill_manual(values = c("same" = "white", "dryer" = "#D73027", "wetter" = "#DAE2F7"), na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories ERA5-Worldclim, 1970-2000", fill = "Worldclim dryer or wetter compared to ERA5")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")


ggplot(filter(diffcatAI, cat.AI.era5 != cat.AI.cmip6)) + geom_raster(aes(x=lon, y = lat,  fill = era5_cmip6))+
    scale_fill_manual(values = c("same" = "white", "dryer" = "#D73027", "wetter" = "#DAE2F7"), na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories ERA5-CMIP6, 1970-2000", fill = "CMIP6 dryer or wetter compared to ERA5")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")


ggplot(filter(diffcatAI, cat.AI.wdc != cat.AI.cmip6)) + geom_raster(aes(x=lon, y = lat,  fill = wdc_cmip6))+
    scale_fill_manual(values = c("same" = "white", "dryer" = "#D73027", "wetter" = "#DAE2F7"), na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories Worldclim-CMIP6, 1970-2000", fill = "CMIP6 dryer or wetter compared to Worldclim")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")


ggplot(filter(diffcatAI, cat.AI.era5 != cat.AI.wdc)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI.wdc))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories ERA5-Worldclim, 1970-2000", fill = "Worldclim dryer or wetter compared to ERA5")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")


ggplot(filter(diffcatAI, cat.AI.era5 != cat.AI.cmip6)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI.cmip6))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories ERA5-CMIP6, 1970-2000", fill = "CMIP6 dryer or wetter compared to ERA5")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")


ggplot(filter(diffcatAI, cat.AI.wdc != cat.AI.cmip6)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI.cmip6))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories Worldclim-CMIP6, 1970-2000", fill = "CMIP6 dryer or wetter compared to Worldclim")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

3.1.1 Compare ETO and pr in WDC, ERA5, CMIP6


bpr <- c(-400, 0, 400, 800, 1200, 1600, 5000)
colvecrdbu <-  brewer.pal(n = 11, name ="RdYlBu")

col.pr <- c("(-400,0]" = colvecrdbu[5], "(0,400]"= colvecrdbu[6],"(400,800]"= colvecrdbu[7], "(800,1.2e+03]" = colvecrdbu[8], "(1.2e+03,1.6e+03]"= colvecrdbu[9], "(1.6e+03,5e+03]" = colvecrdbu[10])

calib$pr.cat <- calib$pr %>% cut(breaks = bpr) 

g.pr.wdc <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual precipitation in Wordclim", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.pr.era5 <- ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
      scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual precipitation in ERA5", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.pr.cmip6 <- ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual precipitation CMIP6", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

l.pr <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top", legend.key.width = unit(2, "cm"))

#)
bpr <- c(-400, 0, 400, 800, 1200, 1600, 5000)
colvecrdbu <-  brewer.pal(n = 11, name ="RdYlBu")

col.et0 <- c("(-400,0]" = colvecrdbu[8], "(0,400]"= colvecrdbu[7],"(400,800]"= colvecrdbu[6], "(800,1.2e+03]" = colvecrdbu[5], "(1.2e+03,1.6e+03]"= colvecrdbu[4], "(1.6e+03,5e+03]" = colvecrdbu[3])

calib$et0.cat <- calib$ET0 %>% cut(breaks = bpr) 

g.et0.wdc <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual evapotranspiration Wordclim", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.et0.era5 <- ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual evapotranspiration ERA5", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.et0.cmip6 <- ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual evapotranspiration CMIP6", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

l.et0 <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    geom_raster(aes(x = lon, y = lat), fill = "white")+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top", legend.key.width = unit(2, "cm"))
ggarrange(
  ggarrange(plotlist = list(g.pr.wdc,g.pr.era5, g.pr.cmip6, l.pr), nrow = 4, ncol = 1),
  ggarrange(plotlist = list(g.et0.wdc,g.et0.era5, g.et0.cmip6, l.et0), nrow = 4, ncol = 1),
  ncol = 2)

3.2 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000 with Antarctica

3.2.1 Map aridity categories for the 3 databases


calib <- rbind(cmip6vA, era5vA, wdcvA)

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
# col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"
calib$cat.AI <- factor(calib$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "ERA5, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "bottom"),

nrows = 4, ncol = 1

)

3.2.2 Proportion aridity category by subregion, continent et database

calib$same.cat <- NA
calib$cat.ref <- NA

for(i in 1:nrow(calib)){
  loni = calib$lon[i]
  lati = calib$lat[i]
  catrefi <- subset(calib, lon == loni & lat == lati & source == "CMIP6")$cat.AI %>% as.character()
  cati <- calib$cat.AI[i] %>% as.character()
  calib$same.cat[i] <- ifelse(cati == catrefi, "y",
                              ifelse(cati!= catrefi, "n", 
                                     NA))
  calib$cat.ref[i] <- catrefi
}

write.table(calib, "figures_article/calib.txt")
# table_calib_c <- calib %>% select(c("Continent","source","cat.AI")) %>%
#   reshape2::melt(id.vars = c("Continent","source")) %>%
#   group_by(Continent,source,value) %>%
#   summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
#   ungroup()
# 
# table_calib_c$value[is.na(table_calib_c$value) == T] <- "NA"
# 
# ggplot(table_calib_c)+
#   geom_col(aes(x = Continent, y = perc, fill = value, position = "fill"))+
#   scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
#   theme_minimal()+
#   theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
#   facet_grid(rows = vars(source))+
#   labs(x = "Continent", y = "%", fill = "Aridity\ncategory", y = "Proportion of aridity category")
# 


table_calib <- calib %>% select(c("Continent","Name","Acronym","source","cat.AI")) %>%
  reshape2::melt(id.vars = c("Continent","Name","Acronym","source")) %>%
  filter(!is.na(value)) %>%
  group_by(Continent,Name,Acronym,source,value) %>%
  summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
  ungroup() %>% filter(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC"))

table_calib$value[is.na(table_calib$value) == T] <- "NA"
table_calib$Continent[which(table_calib$Continent == "EUROPE-AFRICA")] <- "EUROPE"
table_calib$Continent[which(table_calib$Continent == "CENTRAL-AMERICA")] <- "CENTRAL\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "NORTH-AMERICA")] <- "NORTH\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "SOUTH-AMERICA")] <- "SOUTH\nAMERICA"

table_calib$Continent <- factor(table_calib$Continent, levels = c("AFRICA","EUROPE","ASIA","OCEANIA","NORTH\nAMERICA","CENTRAL\nAMERICA","SOUTH\nAMERICA","POLAR","NA"))

table_calib$value <- factor(table_calib$value, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))

ggplot(table_calib)+
  geom_col(aes(x = Acronym, y = perc, fill = value, position = "fill"))+
  scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))+
  theme_minimal(base_size = 12)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
        legend.position = "bottom", 
        strip.text.x = element_text(size = 6))+
  facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
  labs(x = "", y = "%", fill = "Aridity\ncategory")


ggplot(table_calib)+
  geom_col(aes(x = Acronym, y = count, fill = value, position = "fill"))+
  scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))+
  scale_y_continuous(position = "right")+
  theme_minimal(base_size = 12)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
        legend.position = "bottom", 
        strip.text.x = element_text(size = 6))+
  facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
  labs(x = "", y = "Count", fill = "Aridity\ncategory")

3.2.3 Pie chart

library(ggrepel)

pie_calib <- table_calib %>% group_by(source, value) %>% summarise(gridcells = sum(count)) %>% mutate(percent = round(gridcells/sum(gridcells)*100, 1), lab.y = (rev(cumsum(rev(percent)))) - percent*0.5)


ggplot(pie_calib)+geom_bar(aes(x = "", y = percent, fill = value), width = 1, stat = "identity")+
  coord_polar("y", start = 0)+
  scale_fill_manual(values = col.cat, na.translate = F)+
  facet_grid(cols = vars(source))+
  geom_label_repel(aes(x = "", y = lab.y, label = percent), nudge_x = 0.5)+
  labs(fill = "")+
  theme_void()+theme(legend.position = "bottom")

pie_calib <- rbind(cmip6vA, era5vA, wdcvA)

pie_calib$cat.AI <- factor(pie_calib$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold")) 

pie_calib <- pie_calib %>% 
  subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(c("source","cat.AI")) %>%
   table() %>% as.data.frame() %>%
   group_by(source) %>%
    mutate(perc = round(Freq/sum(Freq)*100, 1),
         lab.y = rev(cumsum(rev(Freq))) - Freq*0.5) 

ggplot(pie_calib)+geom_bar(aes(x = "", y = Freq, fill = cat.AI), width = 1, stat = "identity")+
  coord_polar("y", start = 0)+
  scale_fill_manual(values = col.cat, na.translate = F)+
  facet_grid(cols = vars(source))+
  geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
  labs(fill = "")+
  theme_void()+theme(legend.position = "bottom")


# cmip6 <- read.table("cmip6_AI.txt")
# 
# cmip6.ref <- cmip6 %>% filter(period == "1970_2000") %>% select(c("lon","lat","model","period","source","AI","spr","t","sET0")) %>% setNames(c("lon","lat","model","period","source","AI.ref","spr.ref","t.ref","ET0.ref"))
# cmip6.ref$cat.AI <- cmip6.ref$AI.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
# cmip6.ref$cat.AI[which(cmip6.ref$ET0.ref < 400)] <- "Cold"
# 
# write.table(cmip6.ref, "figures_article/cmip6.ref.txt")

cmip6.ref <- read.table("figures_article/cmip6.ref.txt")

wdc <- filter(calib, source == "Worldclim") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.wdc"))
era5 <- filter(calib, source == "ERA5") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.era5"))
cmip6.cal <- filter(calib, source == "CMIP6") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmip6"))

casesm2 <- filter(cmip6.ref, source == "CAS-ESM2") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.casesm2"))
cesm <- filter(cmip6.ref, source == "CESM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cesm"))
cmcc <- filter(cmip6.ref, source == "CMCC") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmcc"))
cmccesm2 <- filter(cmip6.ref, source == "CMCC-ESM2") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmccesm2"))
cnrm <- filter(cmip6.ref, source == "CNRM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cnrm"))
ecearth3 <- filter(cmip6.ref, source == "EC-Earth3") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.ecearth3"))
fgoals <- filter(cmip6.ref, source == "FGOALS") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.fgoals"))
gfdlesm4 <- filter(cmip6.ref, source == "GFDL-ESM4") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.gfdlesm4"))
inm <- filter(cmip6.ref, source == "INM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.inm"))
inmcm5 <- filter(cmip6.ref, source == "INM-CM5") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.inmcm5"))
mpi <- filter(cmip6.ref, source == "MPI") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.mpi"))
mri <- filter(cmip6.ref, source == "MRI") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.mri"))
noresm <- filter(cmip6.ref, source == "NorESM-2-MM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.noresm"))

all <- merge(wdc, era5, by = c("lon","lat")) %>% 
  merge(cmip6.cal, by = c("lon","lat")) %>% 
  merge(casesm2, by = c("lon","lat")) %>%
  merge(cesm, by = c("lon","lat")) %>%
  merge(cmcc, by = c("lon","lat")) %>%
  merge(cmccesm2, by = c("lon","lat")) %>%
  merge(cnrm, by = c("lon","lat")) %>%
  merge(ecearth3, by = c("lon","lat")) %>%
  merge(fgoals, by = c("lon","lat")) %>%
  merge(gfdlesm4, by = c("lon","lat")) %>%
  merge(inm, by = c("lon","lat")) %>%
  merge(inmcm5, by = c("lon","lat")) %>%
  merge(mpi, by = c("lon","lat")) %>%
  merge(mri, by = c("lon","lat")) %>%
  merge(noresm, by = c("lon","lat"))

all$cmip6.wdc <- with(all, ifelse(cat.AI.cmip6 == cat.AI.wdc, "same","diff"))
all$cmip6.era5 <- with(all, ifelse(cat.AI.cmip6 == cat.AI.era5, "same","diff"))
all$wdc.era5 <- with(all, ifelse(cat.AI.wdc == cat.AI.era5, "same","diff"))


all$casesm2.wdc <- with(all, ifelse(cat.AI.casesm2 == cat.AI.wdc, "same","diff"))
all$cesm.wdc <- with(all, ifelse(cat.AI.cesm == cat.AI.wdc, "same","diff"))
all$cmcc.wdc <- with(all, ifelse(cat.AI.cmcc == cat.AI.wdc, "same","diff"))
all$cmccesm2.wdc <- with(all, ifelse(cat.AI.cmccesm2 == cat.AI.wdc, "same","diff"))
all$cnrm.wdc <- with(all, ifelse(cat.AI.cnrm == cat.AI.wdc, "same","diff"))
all$ecearth3.wdc <- with(all, ifelse(cat.AI.ecearth3 == cat.AI.wdc, "same","diff"))
all$fgoals.wdc <- with(all, ifelse(cat.AI.fgoals == cat.AI.wdc, "same","diff"))
all$gfdlesm4.wdc <- with(all, ifelse(cat.AI.gfdlesm4 == cat.AI.wdc, "same","diff"))
all$inm.wdc <- with(all, ifelse(cat.AI.inm == cat.AI.wdc, "same","diff"))
all$inmcm5.wdc <- with(all, ifelse(cat.AI.inmcm5 == cat.AI.wdc, "same","diff"))
all$mpi.wdc <- with(all, ifelse(cat.AI.mpi == cat.AI.wdc, "same","diff"))
all$mri.wdc <- with(all, ifelse(cat.AI.mri == cat.AI.wdc, "same","diff"))
all$noresm.wdc <- with(all, ifelse(cat.AI.noresm == cat.AI.wdc, "same","diff"))

df1 <- table(all$cmip6.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CMIP6-WDC"))
df2 <- table(all$casesm2.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CASESM2-WDC"))
df3 <- table(all$cesm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CESM-WDC"))
df4 <- table(all$cmcc.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CMCC-WDC"))
df5 <- table(all$cmccesm2.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CMCCESM2-WDC"))
df6 <- table(all$cnrm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CNRM-WDC"))
df7 <- table(all$ecearth3.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "ECEARTH3-WDC"))
df8 <- table(all$fgoals.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "FGOALS-WDC"))
df9 <- table(all$gfdlesm4.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "GFDLESM4-WDC"))
df10 <- table(all$inm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "INM-WDC"))
df11 <- table(all$inmcm5.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "INMCM5-WDC"))
df12 <- table(all$mpi.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "MPI-WDC"))
df13 <- table(all$mri.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "MRI-WDC"))
df14 <- table(all$noresm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "NORESM-WDC"))

df0 <- merge(df1,df2, by = "diffornot") %>% merge(df3, by = "diffornot") %>% 
  merge(df4, by = "diffornot") %>% merge(df5, by = "diffornot") %>% merge(df6, by = "diffornot") %>% 
  merge(df7, by = "diffornot") %>% merge(df8, by = "diffornot") %>% merge(df9, by = "diffornot") %>%
  merge(df10, by = "diffornot") %>% merge(df11, by = "diffornot") %>% merge(df12, by = "diffornot") %>%
  merge(df13, by = "diffornot") %>% merge(df14, by = "diffornot")

prop.diff.wdc <- df0[1,]/(df0[1,] + df0[2,])*100      
prop.diff.wdc[1,1] <- "Percent diff"
                                                      
df.prop.wdc <- rbind(df0, prop.diff.wdc)

kable(df.prop.wdc, caption = "Percent of gridcells in which the aridity category differs from Worldclim", digits = 1) %>%
  kable_styling(bootstrap_options = "bordered")
Percent of gridcells in which the aridity category differs from Worldclim
diffornot CMIP6-WDC CASESM2-WDC CESM-WDC CMCC-WDC CMCCESM2-WDC CNRM-WDC ECEARTH3-WDC FGOALS-WDC GFDLESM4-WDC INM-WDC INMCM5-WDC MPI-WDC MRI-WDC NORESM-WDC
diff 3273.0 4275.0 3557.0 4006.0 6503 4594.0 3691 5083.0 3531.0 4323.0 4318.0 4276.0 3349.0 3843.0
same 18418.0 17407.0 18125.0 17676.0 15179 17097.0 18000 16608.0 18160.0 17359.0 17364.0 17415.0 18333.0 17839.0
Percent diff 15.1 19.7 16.4 18.5 30 21.2 17 23.4 16.3 19.9 19.9 19.7 15.4 17.7

all$casesm2.era5 <- with(all, ifelse(cat.AI.casesm2 == cat.AI.era5, "same","diff"))
all$cesm.era5 <- with(all, ifelse(cat.AI.cesm == cat.AI.era5, "same","diff"))
all$cmcc.era5 <- with(all, ifelse(cat.AI.cmcc == cat.AI.era5, "same","diff"))
all$cmccesm2.era5 <- with(all, ifelse(cat.AI.cmccesm2 == cat.AI.era5, "same","diff"))
all$cnrm.era5 <- with(all, ifelse(cat.AI.cnrm == cat.AI.era5, "same","diff"))
all$ecearth3.era5 <- with(all, ifelse(cat.AI.ecearth3 == cat.AI.era5, "same","diff"))
all$fgoals.era5 <- with(all, ifelse(cat.AI.fgoals == cat.AI.era5, "same","diff"))
all$gfdlesm4.era5 <- with(all, ifelse(cat.AI.gfdlesm4 == cat.AI.era5, "same","diff"))
all$inm.era5 <- with(all, ifelse(cat.AI.inm == cat.AI.era5, "same","diff"))
all$inmcm5.era5 <- with(all, ifelse(cat.AI.inmcm5 == cat.AI.era5, "same","diff"))
all$mpi.era5 <- with(all, ifelse(cat.AI.mpi == cat.AI.era5, "same","diff"))
all$mri.era5 <- with(all, ifelse(cat.AI.mri == cat.AI.era5, "same","diff"))
all$noresm.era5 <- with(all, ifelse(cat.AI.noresm == cat.AI.era5, "same","diff"))

df1 <- table(all$cmip6.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CMIP6-ERA5"))
df2 <- table(all$casesm2.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CASESM2-ERA5"))
df3 <- table(all$cesm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CESM-ERA5"))
df4 <- table(all$cmcc.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CMCC-ERA5"))
df5 <- table(all$cmccesm2.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CMCCESM2-ERA5"))
df6 <- table(all$cnrm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CNRM-ERA5"))
df7 <- table(all$ecearth3.era5) %>% as.data.frame() %>% setNames(c("diffornot", "ECEARTH3-ERA5"))
df8 <- table(all$fgoals.era5) %>% as.data.frame() %>% setNames(c("diffornot", "FGOALS-ERA5"))
df9 <- table(all$gfdlesm4.era5) %>% as.data.frame() %>% setNames(c("diffornot", "GFDLESM4-ERA5"))
df10 <- table(all$inm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "INM-ERA5"))
df11 <- table(all$inmcm5.era5) %>% as.data.frame() %>% setNames(c("diffornot", "INMCM5-ERA5"))
df12 <- table(all$mpi.era5) %>% as.data.frame() %>% setNames(c("diffornot", "MPI-ERA5"))
df13 <- table(all$mri.era5) %>% as.data.frame() %>% setNames(c("diffornot", "MRI-ERA5"))
df14 <- table(all$noresm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "NORESM-ERA5"))

df0 <- merge(df1,df2, by = "diffornot") %>% merge(df3, by = "diffornot") %>% 
  merge(df4, by = "diffornot") %>% merge(df5, by = "diffornot") %>% merge(df6, by = "diffornot") %>% 
  merge(df7, by = "diffornot") %>% merge(df8, by = "diffornot") %>% merge(df9, by = "diffornot") %>%
  merge(df10, by = "diffornot") %>% merge(df11, by = "diffornot") %>% merge(df12, by = "diffornot") %>%
  merge(df13, by = "diffornot") %>% merge(df14, by = "diffornot")

prop.diff.era5 <- df0[1,]/(df0[1,] + df0[2,])*100      
prop.diff.era5[1,1] <- "Percent diff"
                                                      
df.prop.era5 <- rbind(df0, prop.diff.era5)   

kable(df.prop.era5, caption = "Percent of gridcells in which the aridity category differs from Worldclim", digits = 1) %>%
  kable_styling(bootstrap_options = "bordered")
Percent of gridcells in which the aridity category differs from Worldclim
diffornot CMIP6-ERA5 CASESM2-ERA5 CESM-ERA5 CMCC-ERA5 CMCCESM2-ERA5 CNRM-ERA5 ECEARTH3-ERA5 FGOALS-ERA5 GFDLESM4-ERA5 INM-ERA5 INMCM5-ERA5 MPI-ERA5 MRI-ERA5 NORESM-ERA5
diff 3197.0 3934.0 3451.0 3819.0 6969.0 3859.0 3629.0 5402.0 3860.0 4189.0 4015.0 3633.0 2910.0 3683
same 18494.0 17748.0 18231.0 17863.0 14713.0 17832.0 18062.0 16289.0 17831.0 17493.0 17667.0 18058.0 18772.0 17999
Percent diff 14.7 18.1 15.9 17.6 32.1 17.8 16.7 24.9 17.8 19.3 18.5 16.7 13.4 17

df <- table(all$wdc.era5) %>% as.data.frame() %>% setNames(c("diffornot","ERA5-WDC"))
percent.diff <- df[1,] / (df[1,]+df[2,]) * 100
percent.diff[1,1] <- "Percent diff"
df.prop.wdc.era5 <- rbind(df, percent.diff)

kable(df.prop.wdc.era5, caption = "Percent of gridcells in which the aridity category differs between Worldclim and ERA5", digits = 1) %>% 
  kable_styling(bootstrap_options = "bordered")
Percent of gridcells in which the aridity category differs between Worldclim and ERA5
diffornot ERA5-WDC
diff 2850.0
same 18841.0
Percent diff 13.1

4 Number of inhabitants by region

https://earthdata.nasa.gov/data/catalog/sedac-ciesin-sedac-gpwv4-popcount-r11-4.11

density <- raster("Population_density/gpw_v4_population_density_rev11_2020_1_deg.tif")
plot(density)

popcount <- raster("Population_density/gpw_v4_population_count_rev11_2020_1_deg.tif")
plot(popcount)


land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land

popcount.df <- popcount %>%
  projectRaster(to = land_mask) %>% 
  as.data.frame(xy = T) %>% setNames(c("lon","lat","pop")) %>% 
  merge(land_mask.df, by = c("lon","lat")) %>% filter(lm == 1)

era5vApop <- merge(era5vA, popcount.df, by = c("lon","lat")) 

tab.pop.by.cat <- era5vApop %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(cat.AI, Continent) %>%
  summarise(pop = sum(pop)) 

write.csv(tab.pop.by.cat, "population.count.era5.csv")
knitr::kable(tab.pop.by.cat)
cat.AI Continent pop
Arid AFRICA 9.407395e+07
Arid ASIA 2.608739e+08
Arid CENTRAL-AMERICA 1.102838e+07
Arid EUROPE 3.523270e+05
Arid EUROPE-AFRICA 3.354203e+07
Arid NORTH-AMERICA 8.275538e+05
Arid OCEANIA 2.483547e+05
Arid SOUTH-AMERICA 2.714742e+06
Cold ASIA 9.005723e+06
Cold EUROPE 6.822776e+06
Cold NORTH-AMERICA 8.555420e+05
Cold POLAR NA
Cold SOUTH-AMERICA 6.222267e+04
Dry subhumid AFRICA 1.131375e+08
Dry subhumid ASIA 3.724373e+08
Dry subhumid CENTRAL-AMERICA 2.847546e+07
Dry subhumid EUROPE 5.195230e+06
Dry subhumid EUROPE-AFRICA 5.933504e+07
Dry subhumid NORTH-AMERICA 9.392606e+06
Dry subhumid OCEANIA 7.027896e+05
Dry subhumid SOUTH-AMERICA 1.297947e+07
Humid AFRICA 6.854010e+08
Humid ASIA 3.055659e+09
Humid CENTRAL-AMERICA 1.462475e+08
Humid EUROPE 4.999673e+08
Humid EUROPE-AFRICA 1.452014e+08
Humid NORTH-AMERICA 2.915502e+08
Humid OCEANIA 1.529713e+07
Humid POLAR 4.329376e+02
Humid SOUTH-AMERICA 3.724919e+08
Hyperarid AFRICA 5.328568e+07
Hyperarid ASIA 1.970440e+07
Hyperarid EUROPE-AFRICA 4.146658e+07
Hyperarid SOUTH-AMERICA 5.088920e+05
Semi-arid AFRICA 2.154724e+08
Semi-arid ASIA 5.713649e+08
Semi-arid CENTRAL-AMERICA 3.973955e+07
Semi-arid EUROPE 1.841667e+06
Semi-arid EUROPE-AFRICA 8.619193e+07
Semi-arid NORTH-AMERICA 2.806716e+07
Semi-arid OCEANIA 3.403064e+06
Semi-arid SOUTH-AMERICA 2.872740e+07
NA AFRICA 3.321493e+04
NA ASIA 1.762974e+07
NA CENTRAL-AMERICA 7.132541e+05
NA EUROPE 1.481967e+05
NA EUROPE-AFRICA 1.876606e+06
NA NORTH-AMERICA 1.403549e+03
NA OCEANIA 3.896168e+04
NA POLAR NA
NA SOUTH-AMERICA 2.214226e+05

5 Comparison of variables datasets

vars <- c("lon","lat","Continent","Name","Acronym","tas","pr","sfcWind","Rn","ea","ET0","AI","cat.AI")
merge.vars <- c("lon","lat","Continent","Name","Acronym")
calib.wide <- select(era5vA, vars) %>%
               merge(select(wdcvA, vars), 
                     by = merge.vars) %>%
          merge(select(cmip6vA, vars), 
                by = merge.vars) %>%
          setNames(c(merge.vars,
                     "tas.era5","pr.era5", "sfcWind.era5","Rn.era5","ea.era5","ET0.era5","AI.era5","cat.AI.era5",
                     "tas.wdc","pr.wdc", "sfcWind.wdc","Rn.wdc","ea.wdc","ET0.wdc","AI.wdc","cat.AI.wdc",
                     "tas.cmip6","pr.cmip6", "sfcWind.cmip6","Rn.cmip6","ea.cmip6","ET0.cmip6","AI.cmip6","cat.AI.cmip6"))

5.1 Scatterplots with Antarctica

5.1.1 ERA5 compared to WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

5.1.2 ERA5 and WDC

ggarrange(

ggplot(calib.wide,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
  scale_color_manual(values = "#A50026")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = "#A50026")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = "#A50026")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = "#A50026")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = "#A50026")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = "#A50026")+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

knitr::kable(data.frame(precipitation = cor(calib.wide$pr.wdc, calib.wide$pr.era5, use = "pairwise.complete.obs"),
                        temperature = cor(calib.wide$tas.wdc, calib.wide$tas.era5, use = "pairwise.complete.obs"),
                        wind = cor(calib.wide$sfcWind.wdc, calib.wide$sfcWind.era5, use = "pairwise.complete.obs"),
                        Rn = cor(calib.wide$Rn.wdc, calib.wide$Rn.era5, use = "pairwise.complete.obs"),
                        ea = cor(calib.wide$ea.wdc, calib.wide$ea.era5, use = "pairwise.complete.obs"),
                        ET0 = cor(calib.wide$ET0.wdc, calib.wide$ET0.era5, use = "pairwise.complete.obs")), 
             digits = 2, caption = "Correlation coefficients between ERA5 and Worldclim") %>%
  kableExtra::kable_styling(bootstrap_options = "bordered")
Correlation coefficients between ERA5 and Worldclim
precipitation temperature wind Rn ea ET0
0.89 0.99 0.68 0.67 1 0.98
ggplot(calib.wide) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.wide) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.1.3 CMIP6 and WDC

ggarrange(

ggplot(calib.wide,aes(x = tas.wdc))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  scale_color_manual(values = "#FDAE61")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = "#FDAE61")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  scale_color_manual(values = "#FDAE61")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
   scale_color_manual(values = "#FDAE61")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
   scale_color_manual(values = "#FDAE61")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
   scale_color_manual(values = "#FDAE61")+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2
)

5.1.4 CMIP6 and ERA5


ggarrange(
  
ggplot(calib.wide,aes(x = tas.era5))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
  scale_color_manual(values = "#4575B4")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.era5))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
  scale_color_manual(values = "#4575B4")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.era5))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
  scale_color_manual(values = "#4575B4")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.era5))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
   stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = "#4575B4")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.era5))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
   scale_color_manual(values = "#4575B4")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.era5))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = "#4575B4")+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2)

5.2 Scatterplots without Antarctica

calib.noant <- calib.wide %>% filter(lat > -55)

5.2.1 Boxplots


library(ggtext)

ggarrange(
  
ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = pr.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = pr.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = pr.cmip6, fill = "CMIP6"))+
  labs(x = "Mean annual precipitation, mm", y = "", fill = "")+
  xlim(0,6000)+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = tas.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = tas.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = tas.cmip6, fill = "CMIP6"))+
  labs(x = "Surface temperature, °C", y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = sfcWind.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = sfcWind.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = sfcWind.cmip6, fill = "CMIP6"))+
  labs(x = expression("Surface wind speed, " * m~s^{-1}), y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = Rn.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = Rn.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = Rn.cmip6, fill = "CMIP6"))+
  labs(x = expression("Rn - G, " * MJ~m^{-2}~day^{-1}), y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = ea.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = ea.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = ea.cmip6, fill = "CMIP6"))+
  labs(x = "ea, kPa", y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = ET0.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = ET0.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = ET0.cmip6, fill = "CMIP6"))+
  labs(x = expression("ET0, " * mm~day^{-1}), y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ncol = 3, nrow = 2, common.legend = T, legend = "bottom"
)


library(SpatialPack)


coords <- calib.noant[,c("lon","lat")]

df.era5.wdc <- with(calib.noant, 
                    data.frame(precipitation = c(cor(pr.wdc, pr.era5, use = "pairwise.complete.obs"), cor.test(pr.wdc, pr.era5, use = "pairwise.complete.obs")[[3]], cor(pr.wdc, pr.era5, use = "pairwise.complete.obs")^2, modified.ttest(pr.wdc, pr.era5, coords)$corr, modified.ttest(pr.wdc, pr.era5, coords)$p.value),
                        temperature = c(cor(tas.wdc, tas.era5, use = "pairwise.complete.obs"), cor.test(tas.wdc, tas.era5, use = "pairwise.complete.obs")[[3]], cor(tas.wdc, tas.era5, use = "pairwise.complete.obs")^2, modified.ttest(tas.wdc, tas.era5, coords)$corr, modified.ttest(tas.wdc, tas.era5, coords)$p.value),
                        wind = c(cor(sfcWind.wdc, sfcWind.era5, use = "pairwise.complete.obs"), cor.test(sfcWind.wdc, sfcWind.era5, use = "pairwise.complete.obs")[[3]], cor(sfcWind.wdc, sfcWind.era5, use = "pairwise.complete.obs")^2, modified.ttest(sfcWind.wdc, sfcWind.era5, coords)$corr, modified.ttest(sfcWind.wdc, sfcWind.era5, coords)$p.value),
                        Rn = c(cor(Rn.wdc, Rn.era5, use = "pairwise.complete.obs"), cor.test(Rn.wdc, Rn.era5, use = "pairwise.complete.obs")[[3]], cor(Rn.wdc, Rn.era5, use = "pairwise.complete.obs")^2, modified.ttest(Rn.wdc, Rn.era5, coords)$corr, modified.ttest(Rn.wdc, Rn.era5, coords)$p.value),
                        ea = c(cor(ea.wdc, ea.era5, use = "pairwise.complete.obs"), cor.test(ea.wdc, ea.era5, use = "pairwise.complete.obs")[[3]], cor(ea.wdc, ea.era5, use = "pairwise.complete.obs")^2, modified.ttest(ea.wdc, ea.era5, coords)$corr, modified.ttest(ea.wdc, ea.era5, coords)$p.value),
                        ET0 = c(cor(ET0.wdc, ET0.era5, use = "pairwise.complete.obs"), cor.test(ET0.wdc, ET0.era5, use = "pairwise.complete.obs")[[3]], cor(ET0.wdc, ET0.era5, use = "pairwise.complete.obs")^2, modified.ttest(ET0.wdc, ET0.era5, coords)$corr, modified.ttest(ET0.wdc, ET0.era5, coords)$p.value)))

df.era5.cmip6 <- with(calib.noant, 
                    data.frame(precipitation = c(cor(pr.cmip6, pr.era5, use = "pairwise.complete.obs"), cor.test(pr.cmip6, pr.era5, use = "pairwise.complete.obs")[[3]], cor(pr.cmip6, pr.era5, use = "pairwise.complete.obs")^2, modified.ttest(pr.cmip6, pr.era5, coords)$corr, modified.ttest(pr.cmip6, pr.era5, coords)$p.value),
                        temperature = c(cor(tas.cmip6, tas.era5, use = "pairwise.complete.obs"), cor.test(tas.cmip6, tas.era5, use = "pairwise.complete.obs")[[3]], cor(tas.cmip6, tas.era5, use = "pairwise.complete.obs")^2, modified.ttest(tas.cmip6, tas.era5, coords)$corr, modified.ttest(tas.cmip6, tas.era5, coords)$p.value),
                        wind = c(cor(sfcWind.cmip6, sfcWind.era5, use = "pairwise.complete.obs"), cor.test(sfcWind.cmip6, sfcWind.era5, use = "pairwise.complete.obs")[[3]], cor(sfcWind.cmip6, sfcWind.era5, use = "pairwise.complete.obs")^2, modified.ttest(sfcWind.cmip6, sfcWind.era5, coords)$corr, modified.ttest(sfcWind.cmip6, sfcWind.era5, coords)$p.value),
                        Rn = c(cor(Rn.cmip6, Rn.era5, use = "pairwise.complete.obs"), cor.test(Rn.cmip6, Rn.era5, use = "pairwise.complete.obs")[[3]], cor(Rn.cmip6, Rn.era5, use = "pairwise.complete.obs")^2, modified.ttest(Rn.cmip6, Rn.era5, coords)$corr, modified.ttest(Rn.cmip6, Rn.era5, coords)$p.value),
                        ea = c(cor(ea.cmip6, ea.era5, use = "pairwise.complete.obs"), cor.test(ea.cmip6, ea.era5, use = "pairwise.complete.obs")[[3]], cor(ea.cmip6, ea.era5, use = "pairwise.complete.obs")^2, modified.ttest(ea.cmip6, ea.era5, coords)$corr, modified.ttest(ea.cmip6, ea.era5, coords)$p.value),
                        ET0 = c(cor(ET0.cmip6, ET0.era5, use = "pairwise.complete.obs"), cor.test(ET0.cmip6, ET0.era5, use = "pairwise.complete.obs")[[3]], cor(ET0.cmip6, ET0.era5, use = "pairwise.complete.obs")^2, modified.ttest(ET0.cmip6, ET0.era5, coords)$corr, modified.ttest(ET0.cmip6, ET0.era5, coords)$p.value)))

df.wdc.cmip6 <- with(calib.noant, 
                    data.frame(precipitation = c(cor(pr.cmip6, pr.wdc, use = "pairwise.complete.obs"), cor.test(pr.cmip6, pr.wdc, use = "pairwise.complete.obs")[[3]], cor(pr.cmip6, pr.wdc, use = "pairwise.complete.obs")^2, modified.ttest(pr.cmip6, pr.wdc, coords)$corr, modified.ttest(pr.cmip6, pr.wdc, coords)$p.value),
                        temperature = c(cor(tas.cmip6, tas.wdc, use = "pairwise.complete.obs"), cor.test(tas.cmip6, tas.wdc, use = "pairwise.complete.obs")[[3]], cor(tas.cmip6, tas.wdc, use = "pairwise.complete.obs")^2, modified.ttest(tas.cmip6, tas.wdc, coords)$corr, modified.ttest(tas.cmip6, tas.wdc, coords)$p.value),
                        wind = c(cor(sfcWind.cmip6, sfcWind.wdc, use = "pairwise.complete.obs"), cor.test(sfcWind.cmip6, sfcWind.wdc, use = "pairwise.complete.obs")[[3]], cor(sfcWind.cmip6, sfcWind.wdc, use = "pairwise.complete.obs")^2, modified.ttest(sfcWind.cmip6, sfcWind.wdc, coords)$corr, modified.ttest(sfcWind.cmip6, sfcWind.wdc, coords)$p.value),
                        Rn = c(cor(Rn.cmip6, Rn.wdc, use = "pairwise.complete.obs"), cor.test(Rn.cmip6, Rn.wdc, use = "pairwise.complete.obs")[[3]], cor(Rn.cmip6, Rn.wdc, use = "pairwise.complete.obs")^2, modified.ttest(Rn.cmip6, Rn.wdc, coords)$corr, modified.ttest(Rn.cmip6, Rn.wdc, coords)$p.value),
                        ea = c(cor(ea.cmip6, ea.wdc, use = "pairwise.complete.obs"), cor.test(ea.cmip6, ea.wdc, use = "pairwise.complete.obs")[[3]], cor(ea.cmip6, ea.wdc, use = "pairwise.complete.obs")^2, modified.ttest(ea.cmip6, ea.wdc, coords)$corr, modified.ttest(ea.cmip6, ea.wdc, coords)$p.value),
                        ET0 = c(cor(ET0.cmip6, ET0.wdc, use = "pairwise.complete.obs"), cor.test(ET0.cmip6, ET0.wdc, use = "pairwise.complete.obs")[[3]], cor(ET0.cmip6, ET0.wdc, use = "pairwise.complete.obs")^2, modified.ttest(ET0.cmip6, ET0.wdc, coords)$corr, modified.ttest(ET0.cmip6, ET0.wdc, coords)$p.value)))

df.combined <- rbind(
  cbind(Metric = c("r","p-value","r2","Dutilleul","p-value"), df.era5.wdc),
  cbind(Metric = c("r","p-value","r2","Dutilleul","p-value"), df.era5.cmip6),
  cbind(Metric = c("r","p-value","r2","Dutilleul","p-value"), df.wdc.cmip6)
)

knitr::kable(df.combined, 
             digits = 2, caption = "Pearson correlation coefficients and r2 between ERA5, Worldclim and CMIP6") %>%
  kableExtra::kable_styling(bootstrap_options = "bordered") %>%
  kableExtra::group_rows(group_label = "ERA5 vs WDC", start_row = 1, end_row = 5) %>%
  kableExtra::group_rows("ERA5 vs CMIP6", 6,10) %>%
  kableExtra::group_rows("WDC vs CMIP6", 11,15)
Pearson correlation coefficients and r2 between ERA5, Worldclim and CMIP6
Metric precipitation temperature wind Rn ea ET0
ERA5 vs WDC
r 0.88 1.00 0.77 0.89 1.00 0.97
p-value 0.00 0.00 0.00 0.00 0.00 0.00
r2 0.77 1.00 0.59 0.79 0.99 0.94
Dutilleul 0.88 1.00 0.77 0.89 1.00 0.97
p-value 0.00 0.00 0.00 0.00 0.00 0.00
ERA5 vs CMIP6
r 0.89 1.00 0.88 0.99 0.99 0.98
p-value 0.00 0.00 0.00 0.00 0.00 0.00
r2 0.79 1.00 0.77 0.97 0.98 0.96
Dutilleul 0.89 1.00 0.88 0.99 0.99 0.98
p-value 0.00 0.00 0.00 0.00 0.00 0.00
WDC vs CMIP6
r 0.89 1.00 0.78 0.91 0.99 0.96
p-value 0.00 0.00 0.00 0.00 0.00 0.00
r2 0.79 0.99 0.61 0.82 0.97 0.91
Dutilleul 0.89 1.00 0.78 0.91 0.99 0.96
p-value 0.00 0.00 0.00 0.00 0.00 0.00

5.2.2 ERA5 compared to WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

ggarrange(

ggplot(calib.noant,aes(x = tas.wdc))+geom_point(aes( y = tas.era5), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5")+
  lims(x = c(-30, 31), y = c(-30,31))+
   annotate(geom = "text", x =  min(c(-30, 31)) + dist(c(-30, 31))*0.15, y =  max(c(-30, 31)) - dist(c(-30, 31))*0.15, 
            label = paste("r = ", round(cor(calib.noant$tas.wdc, calib.noant$tas.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  lims(x = c(-1, 6000), y = c(-1, 6000))+
   annotate(geom = "text", x =  min(c(-1, 6000)) + dist(c(-1, 6000))*0.15, y =  max(c(-1, 6000)) - dist(c(-1, 6000))*0.15, 
            label = paste("r = ", round(cor(calib.noant$pr.wdc, calib.noant$pr.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5")+
  lims(x = c(0,10), y = c(0,10))+
   annotate(geom = "text", x =  min(c(0,10)) + dist(c(0,10))*0.15, y =  max(c(0,10)) - dist(c(0,10))*0.15,
            label = paste("r = ", round(cor(calib.noant$sfcWind.wdc, calib.noant$sfcWind.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5), alpha = 0.1, shape = 16, size = 0.5, col = colorvec[1])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   lims(x = c(-8,20), y = c(-8,20))+
   annotate(geom = "text", x =  min(c(-8,20)) + dist(c(-8,20))*0.15, y =  max(c(-8,20)) - dist(c(-8,20))*0.15,
            label = paste("r = ", round(cor(calib.noant$Rn.wdc, calib.noant$Rn.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC", y = "ea, kPa, ERA5", col = "Dataset")+
   lims(x = c(0,3), y = c(0,3))+
   annotate(geom = "text",  x =  min(c(0,3)) + dist(c(0,3))*0.15, y =  max(c(0,3)) - dist(c(0,3))*0.15,
            label = paste("r = ", round(cor(calib.noant$ea.wdc, calib.noant$ea.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
lims(x = c(-261,4200), y = c(-261,4200))+
   annotate(geom = "text",  x =  min(c(-261,4200)) + dist(c(-261,4200))*0.15, y =  max(c(-261,4200)) - dist(c(-261,4200))*0.15,
            label = paste("r = ", round(cor(calib.noant$ET0.wdc, calib.noant$ET0.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

ggplot(calib.noant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.noant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.2.3 CMIP6 and WDC

ggarrange(

ggplot(calib.noant,aes(x = tas.wdc))+geom_point(aes( y = tas.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  lims(x = c(-30, 31), y = c(-30,31))+
   annotate(geom = "text", x =  min(c(-30, 31)) + dist(c(-30, 31))*0.15, y =  max(c(-30, 31)) - dist(c(-30, 31))*0.15, 
            label = paste("r = ", round(cor(calib.noant$tas.wdc, calib.noant$tas.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4], colorvec[4]))+
  lims(x = c(-1, 6000), y = c(-1, 6000))+
   annotate(geom = "text", x =  min(c(-1, 6000)) + dist(c(-1, 6000))*0.15, y =  max(c(-1, 6000)) - dist(c(-1, 6000))*0.15, 
            label = paste("r = ", round(cor(calib.noant$pr.wdc, calib.noant$pr.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  lims(x = c(0,10), y = c(0,10))+
   annotate(geom = "text", x =  min(c(0,10)) + dist(c(0,10))*0.15, y =  max(c(0,10)) - dist(c(0,10))*0.15,
            label = paste("r = ", round(cor(calib.noant$sfcWind.wdc, calib.noant$sfcWind.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.cmip6), alpha = 0.1, shape = 16, size = 0.5, col = colorvec[4])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6", col = "Dataset")+
   lims(x = c(-8,20), y = c(-8,20))+
   annotate(geom = "text", x =  min(c(-8,20)) + dist(c(-8,20))*0.15, y =  max(c(-8,20)) - dist(c(-8,20))*0.15,
            label = paste("r = ", round(cor(calib.noant$Rn.wdc, calib.noant$Rn.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6", col = "Dataset")+
   lims(x = c(0,3), y = c(0,3))+
   annotate(geom = "text",  x =  min(c(0,3)) + dist(c(0,3))*0.15, y =  max(c(0,3)) - dist(c(0,3))*0.15,
            label = paste("r = ", round(cor(calib.noant$ea.wdc, calib.noant$ea.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6", col = "Dataset")+
lims(x = c(-261,4200), y = c(-261,4200))+
   annotate(geom = "text",  x =  min(c(-261,4200)) + dist(c(-261,4200))*0.15, y =  max(c(-261,4200)) - dist(c(-261,4200))*0.15,
            label = paste("r = ", round(cor(calib.noant$ET0.wdc, calib.noant$ET0.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.2.4 CMIP6 and ERA5


ggarrange(
  
ggplot(calib.noant,aes(x = tas.era5))+geom_point(aes( y = tas.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6")+
  lims(x = c(-30, 31), y = c(-30,31))+
   annotate(geom = "text", x =  min(c(-30, 31)) + dist(c(-30, 31))*0.15, y =  max(c(-30, 31)) - dist(c(-30, 31))*0.15, 
            label = paste("r = ", round(cor(calib.noant$tas.era5, calib.noant$tas.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.era5))+geom_point(aes( y = pr.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4], colorvec[4]))+
  lims(x = c(-1, 6000), y = c(-1, 6000))+
   annotate(geom = "text", x =  min(c(-1, 6000)) + dist(c(-1, 6000))*0.15, y =  max(c(-1, 6000)) - dist(c(-1, 6000))*0.15, 
            label = paste("r = ", round(cor(calib.noant$pr.era5, calib.noant$pr.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.era5))+geom_point(aes( y = sfcWind.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6")+
  lims(x = c(0,10), y = c(0,10))+
   annotate(geom = "text", x =  min(c(0,10)) + dist(c(0,10))*0.15, y =  max(c(0,10)) - dist(c(0,10))*0.15,
            label = paste("r = ", round(cor(calib.noant$sfcWind.era5, calib.noant$sfcWind.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.era5))+geom_point(aes( y = Rn.cmip6), alpha = 0.1, shape = 16, size = 0.5, col = colorvec[8])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, MJ/m2/day, CMIP6", col = "Dataset")+
   lims(x = c(-8,20), y = c(-8,20))+
   annotate(geom = "text", x =  min(c(-8,20)) + dist(c(-8,20))*0.15, y =  max(c(-8,20)) - dist(c(-8,20))*0.15,
            label = paste("r = ", round(cor(calib.noant$Rn.era5, calib.noant$Rn.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.era5))+geom_point(aes(y = ea.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5, °C", y = "ea, kPa, CMIP6", col = "Dataset")+
   lims(x = c(0,3), y = c(0,3))+
   annotate(geom = "text",  x =  min(c(0,3)) + dist(c(0,3))*0.15, y =  max(c(0,3)) - dist(c(0,3))*0.15,
            label = paste("r = ", round(cor(calib.noant$ea.era5, calib.noant$ea.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.era5))+geom_point(aes(y = ET0.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
lims(x = c(-261,4200), y = c(-261,4200))+
   annotate(geom = "text",  x =  min(c(-261,4200)) + dist(c(-261,4200))*0.15, y =  max(c(-261,4200)) - dist(c(-261,4200))*0.15,
            label = paste("r = ", round(cor(calib.noant$ET0.era5, calib.noant$ET0.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2

)

5.3 Scatterplots only Antarctica

calib.ant <- calib.wide %>% filter(lat < -60)

5.3.1 ERA5 and WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

ggarrange(

ggplot(calib.ant,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

ggplot(calib.ant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.ant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.3.2 CMIP6 and WDC

ggarrange(

ggplot(calib.ant,aes(x = tas.wdc))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2
)

5.3.3 CMIP6 and ERA5


ggarrange(
  
ggplot(calib.ant,aes(x = tas.era5))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
  scale_color_manual(values = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.era5))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.era5))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.era5))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
   stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.era5))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.era5))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2)

5.4 Aridity category difference

5.4.1 ERA5 and WDC

calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.era5 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.era5, cat.AI.wdc, sep = " to "))
table(calib.wide$change.cat.wdc)
## 
##                 Arid to Arid                Arid to Humid 
##                         1359                            1 
##            Arid to Hyperarid            Arid to Semi-arid 
##                           18                          155 
##                 Cold to Arid                 Cold to Cold 
##                           19                        10101 
##         Cold to Dry subhumid                Cold to Humid 
##                            9                          290 
##            Cold to Semi-arid         Dry subhumid to Arid 
##                           20                           23 
##         Dry subhumid to Cold Dry subhumid to Dry subhumid 
##                            5                          205 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          209                          271 
##                Humid to Arid                Humid to Cold 
##                           30                          520 
##        Humid to Dry subhumid               Humid to Humid 
##                          347                         5436 
##           Humid to Semi-arid            Hyperarid to Arid 
##                          210                          230 
##            Hyperarid to Cold       Hyperarid to Hyperarid 
##                            2                          726 
##                   NA to Cold            Semi-arid to Arid 
##                          570                          181 
##            Semi-arid to Cold    Semi-arid to Dry subhumid 
##                            3                          173 
##           Semi-arid to Humid       Semi-arid to Semi-arid 
##                           71                         1077
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.wdc == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Arid to Cold", "dryer",       
                                ifelse(change.cat.wdc == "Cold to Humid", "dryer",
                                ifelse(change.cat.wdc == "Cold to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Cold","dryer",
                                ifelse(change.cat.wdc == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Humid", "dryer",
                                NA))))))))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "ERA5 compared to WDC", title = "Change of category")+
  ylim(-55,90)+
  theme_void(base_size = 15)+
  theme(legend.position = "bottom"),

ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.era5-AI.wdc))+
  binned_scale(aesthetics = "fill", breaks = c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), 
               palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")),
               guide = guide_legend(label.theme = element_text(angle = 0))
               )+
  labs(title = "Changes in AI", fill = "ERA5 compared to WDC")+
  theme_void(base_size = 15)+ylim(-55,90)+
  theme(legend.position = "bottom"),
ncol = 2)

ggplot() + 
  geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc < 277 & pr.era5 - pr.wdc > - 98), aes(x=lon, y = lat,  fill = pr.era5-pr.wdc))+
  geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc > 277), aes(x=lon, y = lat), fill = "#44FF44")+
   geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc < - 98), aes(x=lon, y = lat), fill = "#FF4444")+
  labs(title = "Most extreme differences in precipitation", subtitle = "In red, the 10% most extreme gricells in which Worldclim is wetter than ERA5; \nin green, the 10% most extreme gridcells in which ERA5 is wetter than Worldclim", fill = "precipitation in ERA5 - precipitations in WDC,\nmm/year")+
  theme_void(base_size = 10)+
  theme(legend.position = "bottom")

5.4.2 CMIP6 and WDC

calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.cmip6, cat.AI.wdc, sep = " to "))
table(calib.wide$change.cat.wdc)
## 
##                 Arid to Arid         Arid to Dry subhumid 
##                         1099                            9 
##            Arid to Hyperarid            Arid to Semi-arid 
##                           83                          150 
##                 Cold to Arid                 Cold to Cold 
##                           25                        10491 
##         Cold to Dry subhumid                Cold to Humid 
##                            5                          588 
##            Cold to Semi-arid         Dry subhumid to Arid 
##                           23                           36 
## Dry subhumid to Dry subhumid        Dry subhumid to Humid 
##                          147                          196 
##    Dry subhumid to Hyperarid    Dry subhumid to Semi-arid 
##                            1                          350 
##                Humid to Arid                Humid to Cold 
##                           77                          138 
##        Humid to Dry subhumid               Humid to Humid 
##                          427                         5112 
##           Humid to Hyperarid           Humid to Semi-arid 
##                            1                          435 
##            Hyperarid to Arid            Hyperarid to Cold 
##                           85                            2 
##       Hyperarid to Hyperarid                   NA to Cold 
##                          652                          570 
##            Semi-arid to Arid    Semi-arid to Dry subhumid 
##                          520                          146 
##           Semi-arid to Humid       Semi-arid to Hyperarid 
##                          111                            7 
##       Semi-arid to Semi-arid 
##                          775
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.wdc == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Arid to Cold", "dryer",       
                                ifelse(change.cat.wdc == "Cold to Humid", "dryer",
                                ifelse(change.cat.wdc == "Cold to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Cold","dryer",
                                ifelse(change.cat.wdc == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Humid", "dryer",
                                NA))))))))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "CMIP6 compared to WDC", "Change of category")+
  ylim(-55,90)+
  theme_void(base_size = 15)+
  theme(legend.position = "bottom"),

ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.cmip6-AI.wdc))+
   borders()+
  binned_scale(aesthetics = "fill", breaks = c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")), 
                                                                                                              guide = guide_legend(label.theme = element_text(angle = 0)))+
  labs(title = "Change in AI", fill = "CMIP6 compared to WDC")+
  theme_void(base_size = 15)+ylim(-55,90)+
  theme(legend.position = "bottom"),

ncol = 2)

5.4.3 CMIP6 and ERA5

calib.wide$same.cat.era5 <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.era5, "y","n"))
calib.wide$change.cat.era5 <- with(calib.wide, paste(cat.AI.cmip6, cat.AI.era5, sep = " to "))
table(calib.wide$change.cat.era5)
## 
##                 Arid to Arid         Arid to Dry subhumid 
##                          989                            6 
##                Arid to Humid            Arid to Hyperarid 
##                            2                          230 
##            Arid to Semi-arid                 Cold to Cold 
##                          114                        10377 
##         Cold to Dry subhumid                Cold to Humid 
##                           10                          745 
##         Dry subhumid to Arid Dry subhumid to Dry subhumid 
##                           14                          169 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          219                          328 
##                Humid to Arid                Humid to Cold 
##                            5                           62 
##        Humid to Dry subhumid               Humid to Humid 
##                          382                         5456 
##           Humid to Semi-arid            Hyperarid to Arid 
##                          285                           14 
##       Hyperarid to Hyperarid                     NA to NA 
##                          725                          570 
##            Semi-arid to Arid    Semi-arid to Dry subhumid 
##                          511                          146 
##           Semi-arid to Humid       Semi-arid to Hyperarid 
##                          121                            3 
##       Semi-arid to Semi-arid 
##                          778
calib.wide$dryer.era5 <- with(calib.wide, ifelse(change.cat.era5 == "Arid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.era5 == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.era5 == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.era5 == "Cold to Humid", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.era5 == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to Cold","dryer",
                                ifelse(change.cat.era5 == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.era5 == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.era5 == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.era5 == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.era5 == "Hyperarid to Semi-arid", "dryer",
                                NA)))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.era5 == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.era5))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "CMIP6 compared to ERA5", "Change of aridity category")+
  ylim(-55,90)+
  theme_void()+
  theme(legend.position = "bottom"),


ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.cmip6-AI.era5))+
   borders()+
  binned_scale(aesthetics = "fill", breaks =c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")))+
  labs(title = "Change of AI", fill = "CMIP6 compared to ERA5")+
  theme_void()+ylim(-55,90)+
  theme(legend.position = "bottom"),
ncol = 2)

6 Aridity Index Hargreaves


calib$

7 NDVI

# Download Gimms NDVI
library(gimms)
library(gganimate)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4")
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land

#ndvi.ref <- downloadGimms(x= 1981,y= 2000,dsn = "NDVI")

list.ndvi <- list.files(path = "NDVI", pattern = "ndvi.*\\.nc4$", full.names = T)[1:2]

ndvi.rast <- rasterizeGimms(list.ndvi, keep = 0) %>% mean(na.rm = T) %>% projectRaster(land_mask) # rasterize Gimms does the quality control and removes NAs
ndvi.df <- ndvi.rast %>% as.data.frame(xy = T) %>% setNames(c("lon","lat","ndvi")) %>% merge(land_mask.df, by = c("lon","lat")) %>% filter(lm != 0)

write.table(ndvi.df, "NDVI/ndvi.df.txt")
ndvi.df <- read.table("NDVI/ndvi.df.txt") 

ggplot(subset(ndvi.df, lm == 1))+geom_raster(aes(x=lon, y = lat, fill = ndvi))+
  scale_fill_continuous_sequential(palette = "Greens 3", rev = T)+
  ylim(-55,90)+labs(fill = "NDVI")+
  theme_void()+theme(legend.position = "bottom")


calibn <- merge(calib, ndvi.df, by = c("lon","lat"))

ggplot(subset(calibn, !is.na(cat.AI)))+geom_boxplot(aes(x=cat.AI, y = ndvi, col = cat.AI))+facet_grid(rows = vars(source))+
  scale_color_manual(values = col.cat)+
  labs(x = "", y = "NDVI")+
  theme_minimal()+theme(legend.position = "none")


aov(ndvi~cat.AI, data = group_by(calibn, source)) %>% summary()
##                Df Sum Sq Mean Sq F value Pr(>F)    
## cat.AI          5   1026  205.17    8708 <2e-16 ***
## Residuals   43200   1018    0.02                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 23577 observations effacées parce que manquantes

range.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>% 
  summarise(range.ndvi = paste(round(quantile(ndvi, na.rm = T)[c(2,4)],2), collapse = " to ")) %>%
  cast(source~cat.AI)

cor.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>% 
  summarise(pearson.r = round(cor(ndvi, AI, use = "pairwise.complete"),2)) %>%
  cast(source~cat.AI)

knitr::kable(range.ndvi, caption = "NDVI range by source and aridity category")
NDVI range by source and aridity category
source Cold Hyperarid Arid Semi-arid Dry subhumid Humid
CMIP6 0.3 to 0.47 0.09 to 0.11 0.1 to 0.19 0.17 to 0.36 0.24 to 0.49 0.4 to 0.69
ERA5 0.27 to 0.47 0.09 to 0.11 0.11 to 0.21 0.19 to 0.37 0.24 to 0.48 0.42 to 0.68
Worldclim 0.29 to 0.47 0.09 to 0.11 0.1 to 0.19 0.18 to 0.34 0.28 to 0.44 0.45 to 0.7
knitr::kable(cor.ndvi, caption = "Correlation between NDVI and AI by source and aridity category")
Correlation between NDVI and AI by source and aridity category
source Cold Hyperarid Arid Semi-arid Dry subhumid Humid
CMIP6 -0.08 0.15 0.36 0.38 -0.01 0.31
ERA5 -0.16 0.08 0.47 0.30 0.08 0.38
Worldclim 0.05 0.05 0.49 0.39 0.14 0.52
knitr::knit_exit()

7.0.1 Maps CMIP6 ERA5 // WDC

In red: when CMIP6 or ERA5 are warmer than Wordlclim

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.era5-tas.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: ERA5 - WDC")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - WDC")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when CMIP6 or ERA5 are wetter than Worldclim

cowplot::plot_grid(
  
  ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.era5-pr.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: ERA5 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when wind speed is higher for CMIP6 or ERA5 compared to Worldclim

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.era5-sfcWind.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: ERA5 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

7.0.2 Maps CMIP6 WDC // ERA5

In red: when CMIP6 or WDC are warmer than ERA5

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.wdc-tas.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: WDC - ERA5")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - ERA5")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when CMIP6 or WDC are wetter than ERA5

cowplot::plot_grid(
  
  ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.wdc-pr.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: WDC - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when wind speed is higher for CMIP6 or WDC compared to ERA5

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.wdc-sfcWind.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: WDC - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)